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INTRODUCTION 

In the last few years, major advances have been made in the numerical simula- 
tion of wall-bounded transitional and turbulent shear flows. So far, most of the 
emphasis has been on flows with planar boundaries. Transitional flows in a flat- 
plate boundary layer (ref. 1) and in a channel (ref. 2), as well as turbulent channel 
flow (ref. 3), have been investigated without modeling the near-wall region; Moin 
(ref. 4) gives a critique of these investigations. Except for the study of Patera 
and Orszag (ref. 5) of axisymmetric pipe flow and the simulation of turbulent flow in 
annuli by Schumann (ref. 6), using modeled boundary conditions, very little work has 
been reported on nonplanar flows. 

In this paper we present a new numerical technique for simulating three- 
dimensional, unsteady, incompressible pipe flows and demonstrate its utility and 
accuracy. Each vector function in the expansion of the velocity field is divergence- 
free and satisfies the boundary conditions for viscous flow. Some of the benefits of 
the expansion technique are as follows: (1) pressure is eliminated from the dynam- 

ics, (2) only two unknowns per "mesh point" are required, (3) it provides implicit 
treatment of the viscous terms at no extra computational cost, and (4) no fractional 
time-steps are required. 

In addition, the method uses spectral expansions: Fourier series in the azi- 

muthal (6) and streamwise (x) directions and global polynomials in the radial (r) 
direction. Thus, for smooth velocity fields, we expect rapid convergence of our 
expansions, independent of boundary constraints, as long as the radial polynomials 
are eigenfunctions of a singular Sturm-Liouville problem (ref. 7). In general, 
Chebychev or Legendre polynomials are good candidates, but we show that for cylindri- 
cal geometry a certain choice of the Jacobi polynomials is particularly advantageous 
in minimizing the coupling of the resulting equations for the expansion coefficients 
while satisfying the analytical behavior of the flow variables near r = 0. As dis- 
cussed below, the method has been tested on the linear stability problem for 
Poiseuille flow. 
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MATHEMATICAL FOUNDATIONS 


The governing equations are the incompressible Navier-Stokes equations for the 
velocity _u and the disturbance pressure p. 


1 

^ X _u = -V(p + - u^) 


dP ^ ^ 1 

dx X Re 


V^u 


( 1 ) 


V • u = 0 . (2) 

Here ^ = Vx^ is the vorticity. Re is the Reynolds number, is the constant 
mean pressure gradient, and the density p = 1 everywhere. The boundary condition 
at the pipe wall (r = 1) is assume periodic boundary conditions in 

the X direction with period L. 


We proceed formally by defining the projection operator ^ which projects an 
arbitrary vector field into the space of divergence-free fields satisfying tangency 
at the boundary; that is , if Q is an arbitrary vector field, then we have the 
unique decomposition, 


Q = 


(3) 


where satisfies 


V • i) = 0 


i!i • ill “0 

boundary 


(4) 


and ^ is the operator that accomplishes this decomposition. 




(5) 


By applying this operator to the momentum equation, we obtain the time-evolution of 
the divergence-free velocity field as (ref. 8) 


du 

— ^(V X V X u) = -^(o) X u) - — e 
dt Re — — — dx X 


( 6 ) 


Our overall strategy is to expand jj in terms of divergence-free vector func- 
tions satisfying the viscous boundary conditions and periodicity in x and 0. We 
then substitute this expansion into equation (1) and apply a weighted residual 
method which mimics the application of the projection operator to obtain evolution 
equations for the expansion coefficients. 
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EXPANSION METHOD 


We write the velocity field u as the expansion 

u(r,0,x,t) = exp(ikx + iJie) 

n,k,£ 


( 7 ) 


where each expansion vector satisfies 

7 • [i,(r) exp(ikx + i£6)J = 0 


( 8 ) 


and 


X„(l) = 0 (9) 

We derive a system of ordinary differential equations (ODE's) for the coefficients 
a , by substituting the above expansion into the momentum equation and taking 
the inner product of the result with a set of weight vectors which are divergence- 
free, 

V • [^(r) exp(-ikx - iJi0)J = 0 (10) 


and satisfy the inviscid boundary condition, 


— m 


e = 0 
r 


( 11 ) 


[The (k,^,) dependence of ^ and ^ is suppressed.] This weighted residual method 
mimics the application of the projection operator. For example if 
£ = C (r) exp(-ikx - i£0) then 



Vv> dV = 



O dV + 



* IL) dS 


S 


( 12 ) 


= 0 

for an arbitrary scalar field ip with periodic boundary conditions in x. 
The result for each wave vector (k,£) is a system of ODE's 


where 


A a + 


_L_ 

Re 


B a 


f 


A 


mn 



• Y rdr 


(13) 


(14a) 
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(14b) 


i 


“ fC'VxVxx rdr 
mn / “m 


i 

+ 2ttL £ rdr 


(14c) 


and where denotes double-Fourier transformation in x and 6. Thus, except for 

the nonlinear term, the coupling of the equations occurs only through the radial 
modes . 

The choice of radial functions in ^ and 5 must be made carefully to 
(1) minimize the coupling between radial modes (that is, obtain a banded structure 
for A and B if possible); (2) allow construction of A and B with relative ease; 
(3) obtain effxcxent computation of and (4) obtain rapid convergence while satis— 
the constraints imposed on and including the correct behavior as 

r -»■ 0. 

We find that the sequence of expansion vectors ix • • -I defined 

(“-1 0 0 n n I 

in the following satisfies the above requirements. For n> 0 the (r,6,x) compo- 
nents are 



ikq 


£±1 


kq 


£±1 

n 


_ 1 . £±1 . 
r dr ^n 


£ £±li 

7 'In 


(k / 0) 


(15) 


where 


,0^) 


qj = rl^l(l - r 2 )%^UI)(r 2 ) 


§n shifted Jacobi polynomial (ref. 9), 


(16) 


= P^°’^^2y - 1) 


(17) 


satisfying the orthogonality condition 


1 

/ 


y^ g^*"^ (y) g^^^ (y) dy = 6 

m n m m,n 


(18) 
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For the case k = 0, the above expansion vectors clearly are not complete and 
must be replaced by an alternative set. A convenient choice is given by (n > 0) 



It is a simple matter to verify that the expansion vectors defined above yield 
the correct behavior of ^(r,k,Jl) as r ^ 0; for example, if Z > 0, 


. Z-1 

ug - lyr 

u ^ Br^ 

X 

where y and 3 are complex constants. The additional vector x is required 

because u (k 0) and u (k = 0) would otherwise have a double zero at the’ wall. 

D X 

The corresponding weight vectors are essentially the curl of the x’s. More 
specifically, if k ^ 0 the weight vectors may be expressed as 



while the 


—n 


have the form 


+ 





( 20 ) 


( 21 ) 


As a result, the (+) vectors are uncoupled from the (-) vectors. The resulting 
matrices A and B are nonadiagonal, except for an additional nonzero row and 
a column owing to the vector The limited bandwidth of the viscous 

matrix B results from the particular choice of the polynomials g^^^ given by 
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equation (17). In particular, the Laplacian operator in the (r,0) plane is equiv- 
alent to a tridiagonal matrix in the following sense: 


[q^(r) exp(i£e] = (7 ^ exp(ijL0) 


CONVERGENCE TESTS 

As a test, the method described above was applied to the problem of determining 
the time eigenvalues for linearized flow in a pipe. The calculations were performed 
on a CDC 7600 computer. We assume u ^ exp(Xt) and order the eigenvalues such that 
Real (X^) > Real (X 2 ) > .••• The results for X^ with Re = 9600, £ = 1, k = 1, 
are given in the table below, where N + 2 is the number of radial modes in the 
expansion. 

N X 

20 -0.02312 - i 0.95050 

25 -0.02317074 - i 0.95048142 

30 -0.023170795769 - i 0.950481396659 

35 -0.023170795764 - i 0.950481396668 

Note that the convergence is exponential in N or some power of N, typical of 
spectral methods (ref. 7), and that there is no indication of significant round-off 
errors. The results agree with that of Salwen et al. (ref. 10), who obtained 
Xi = 0.02317 - i(0. 95048), using an expansion in Stokes’ eigenfunctions to solve the 
linear stability problem. 

In figure 1 we show convergence of some of the higher eigenvalues for the case 
k= 1, £= 1, Re= 3000, a Reynolds number in the range where a number of interest- 
ing transitional phenomena have been observed experimentally. Note that a large 
number of eigenvalues are predicted accurately for 30 to 35 radial expansion func- 
tions. In figure 2, the amplitudes of the coefficients of for the eigenvalue 

\l are shown as a function of n for four Reynolds numbers with k = 1 and £ = 1. 
Each of the x normalized so that 


1 



0 
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.iiul I lu' iMclial expansion N = 35 (37 expansion functions) was used for all cases. 
A)*.:! in, Llio coefficients approach zero exponentially in some power of n. 


SUMMARY 

A new numerical method has been developed to investigate three-dimensional, 
unsteady pipe flows using a new velocity-vector expansion method. Each vector func- 
tion in the expansion set is divergence-free and satisfies the boundary conditions 
for viscous flow. Other features of the general technique are as follows: (1) pres- 

sure is eliminated from the dynamics; (2) only two unknox^s per ”mesh point” are 
required; (3) there is rapid convergence of spectral methods; (4) there is implicit 
treatment of the viscous term at no extra computational cost; and (5) no fractional 
time-steps are required. In the present application of the method to flow in a pipe, 
the behavior of each flow variable near the computational singular point is treated 
rigorously and expansions in Jacobi polynomials have been shown to be particularly 
advantageous. The method has been tested on the linear stability problem for 
Poiseuille flow and has demonstrated rapid convergence of the eigenvalues and eigen- 
functions as the number of radial modes is increased. 
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Fig. 1. Convergence of Real (A^); k = J!, = 1, Re = 3000. 



2. Amplitudes of the expansion coefficients for the eigenvector corresponding 
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